# Define scientific model and create simulated data

## Define function for calculating variance of a normal distirbution based on min and max value and set parameters
```{r}
library(rethinking)

# Load dataset from module_1
load("./data/d_sol_cw_w_in.RData")

## Set parameters
a_C <- 2.4
b_C_S <- 0.3
b_C_T <- 0.6
a_T <- 278
b_T_S <- 0.6
```

# 1. Create dataset for assumption m1: T_S <- Sol -> C_MW
```{r}
d_sim_m1 <- d_sol_cw

d_sim_m1$Time <- d_sol_cw$Time

d_sim_m1$Sol <- d_sol_cw$Sol
## Standardize predictor variables
d_sim_m1$Sol_std <- standardize(d_sim_m1$Sol)

d_sim_m1$T_S <- d_sol_cw$T_S
## Standardize predictor variables
d_sim_m1$T_S_std <- standardize(d_sim_m1$T_S)

## Generate simulated C_MW variable
d_sim_m1$C_MW <- rnorm(length(d_sim_m1$Time), a_C + b_C_S*d_sim_m1$Sol_std, 0.4)
d_sim_m1$C_MW_std <- standardize(d_sim_m1$C_MW)

```

# 2. Create dataset for assumption m2: Sol -> T_S -> C_MW
```{r}
d_sim_m2 <- d_sim_m1

## Generate simulated C_MW variable
d_sim_m2$C_MW <- rnorm(length(d_sim_m1$Time), a_C + b_C_T*d_sim_m2$T_S_std, 0.4)
d_sim_m2$C_MW_std <- standardize(d_sim_m2$C_MW)
```

# 3. Create dataset for assumption m3: T_S <- Sol -> C_MW with Sol -> T_S
```{r}
## Copy use Sol values from before
d_sim_m3 <- d_sim_m2

## Generate simulated C_MW variable
d_sim_m3$C_MW <- rnorm(length(d_sim_m1$Time), a_C+b_C_S*d_sim_m3$Sol_std+b_C_T*d_sim_m3$T_S_std, 0.4)
d_sim_m3$C_MW_std <- standardize(d_sim_m3$C_MW)
```

## Generate a noisy dataset

# Plot the simulated radiation
```{r}

windows(); plot(d_sol_cw$C_MW)
windows(); plot(d_sim_m1$C_MW)
windows(); plot(d_sim_m2$C_MW)
windows(); plot(d_sim_m3$C_MW)

windows(); plot(d_sol_cw$T_S_std)
windows(); plot(d_sim_m2$T_S_std)
windows(); plot(d_sim_m1$T_S_std)

windows(); plot(d_sol_cw$Sol_std)
windows(); plot(d_sim_m1$Sol_std)

```

# Export the data for processing in Python
```{r}
save(d_sim_m1, d_sim_m2, d_sim_m3, file = "./data/d_sim_in.RData")
```

## Run all above
```{r}

```